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I.  INTRODUCTION 


The  study  of  time-dependent  flows  past  bluff  bodies  has  historically  been  the  focus 
of  a  great  deal  of  scientific  attention  owing  to  its  relevance  to  many  and  diverse  applications, 
ranging  from  submerged  structures  to  hot  wire  anemometers.  Some  forty  years  ago, 
Keulegan  and  Carpenter  commenced  the  methodical  investigation  of  the  fluid-structure 
interactions  which  occur  when  bodies  are  immersed  in  unsteadily  flowing  fluids.  Today,  the 
effect  of  the  parameters  relevant  to  the  problem,  such  as  Reynolds  number,  Keulegan- 
Carpenter  number  and  relative  roughness,  to  name  a  few,  is  much  better  understood  thanVc 
to  ongoing  research  in  this  area.  In  spite  of  this  progress,  real  time  prediction  of  the  forces 
caused  by  unsteady  flows  on  submerged  objects,  such  as  those  acting  on  an  underwater 
robotic  arm  or  an  offshore  oil  rig,  seriously  challenges  our  current  theoretical  and 
computational  capabilities. 

One  of  the  principal  empirical  tools  used  by  the  engineer  to  solve  the  problems 
described  above,  Morison's  equation,  is  only  reliable  in  predicting  forces  and  moments  in 
highly  idealized  conditions  with  very  low  or  very  high  flow  oscillation  periods. 
Additionally,  most  real-life  problems  involve  flows  in  the  turbulent  regime,  further 
increasing  the  level  of  difficulty  of  both  analytical  and  numerical  solutions. 

Given  the  cost  and  complexity  of  the  laboratory  apparatus  required  to  reproduce  real 
life  flow  conditions,  theoretical  and  experimental  advances  have  been  paralleled  by  a 
considerable  research  effort  in  the  area  of  computational  fluid  dynamics  to  assist  in  the 
solution  of  problems  related  to  bodies  immersed  in  unsteadily  flowing  fluids.  Progress  in 
numerical  techniques  and  the  ever  increasing  power  of  present  day  computers  are  finally 
making  it  possible  to  overcome  the  barrier  historically  imposed  by  the  physical  and 
numerical  instabilities  which  have  caused  the  modelling  of  turbulent  flows  particularly 
challenging  in  the  past. 

Whereas,  until  recently,  numerical  experimentation  has  only  been  successful  in 
determining  flow  patterns,  Strouhal  numbers  and  force  coefficients  for  flow  regimes  within 
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selected  ranges  of  relatively  low  Reynolds  numbers,  it  is  now  possible  to  obtain  useful  data 
for  more  turbulent  flows  exhibiting  extensive  separation.  Furthermore,  commercial  software 
is  becoming  available  which  affords  researchers  much  greater  flexibility  than  the  ad  hoc 
codes  previously  generated  to  solve  very  specific  categories  of  problems,  allowing  for  the 
analysis  of  a  much  broader  gamut  of  flow  conditions. 

Unfortunately,  in  spite  of  such  great  flexibility,  modem  software  does  not  yet  absolve 
the  user  from  having  to  fine-tune  the  code  by  the  judicious  selection  of  parameters, 
numerical  techniques  and  turbulence  models.  Even  then,  experimentally  obtained  results 
cannot  be  exactly  duplicated  in  all  cases. 

If  complete  and  accurate  solutions  are  not  yet  always  achievable,  however, 
approximate  solutions  can  certainly  aid  in  expanding  our  current  level  of  understanding  and 
provide  the  engineer  with  a  valuable  tool  to  improve  design  optimization.  From  the  above 
discussion,  it  can  be  inferred  that  a  benchmark  body  of  reliable  calculations,  performed  using 
the  most  appropriate  set  of  computation  tools,  is  required  to  calibrate  our  current  empirical 
equations  and  experimental  results. 

The  objective  of  this  investigation  is  to  lay  the  foundations  for  such  an  endeavor  and 
improve  our  current  ability  to  predict  the  effects  of  time-dependent  turbulent  flows  over 
circular  cylinders  through  the  use  of  a  state  of  the  art  commercial  code  generated  by  CFD 
Research  Corporation  (CFDRC).  Whenever  possible,  the  results  achieved  have  been 
compared  to  those  obtained  experimentally.  Besides  their  intrinsic  value,  these  results  will 
aid  in  pointing  out  some  of  the  strengths  and  weaknesses  of  the  code  and  hopefully  add  to 
our  understanding  of  the  physics  of  flow  types  not  previously  analyzed. 
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n.  BACKGROUND  STUDIES 


About  twenty  years  after  Keulegan  and  Carpenter  published  their  ground  breaking 
work,  a  comprehensive  series  of  experiments  concerning  sinusoidally  oscillating  flow  about 
smooth  and  rough  bodies  was  performed  at  the  Naval  Postgraduate  School  (Sarpkaya,  1 976). 
This  endeavor  produced  a  wealth  of  practical  results  and  introduced  the  parameter  p  (= 
Re/K)  to  assess  the  influence  of  scale  in  periodic  flows.  Among  its  findings,  the  dependence 
of  the  force-transfer  coefficients  on  Re,  K,  and  k/D  was  shown  graphically  and  amply 
discussed.  As  a  continuation  to  this  research  effort,  cases  were  analyzed  which  involved  a 
coexisting  flow  as  specified  by  U  =  U0  -  Um  cos  (2*t/T),  in  which  U0  represents  the  steady 
mean  velocity  and  Um  the  amplitude  of  sinusoidal  oscillations  (Storm  1984).  Since  the  mid¬ 
eighties,  several  numerical  predictions  of  the  Strouhal  number,  the  pressure  distribution,  and 
the  evolution  of  the  lift  and  drag  forces  in  steady  and  time-dependent  ambient  flows  have 
been  performed  in  an  attempt  to  obtain  results  consistent  with  experimental  data.  Here,  only 
the  more  recent  and  relevant  investigations  will  be  briefly  reviewed. 

A  finite-difference  analysis  of  the  Navier-Stokes  equations  for  a  sinusoidally 
oscillating  ambient  flow  about  a  circular  cylinder  at  K  =  5  (Re  =  1000)  and  K  =  7  (Re  =  700) 
has  been  attempted  by  Baba  &  Miyata  (1987),  assuming  a  physically  unrealistic  symmetric 
wake  in  both  simulations.  The  results  have  shown  that  the  calculations  could  be  carried  out 
only  for  short  times  (less  than  two  cycles  of  flow  oscillation)  with  a  nonsuper  computer. 

A  similar  method  was  used  to  analyze  three  cases  (K  =  5,  7,  and  10)  at  Reynolds 
numbers  around  104  (Murashige  et  al.  1989).  The  flow  was  perturbed  by  artificial  means  to 
trigger  asymmetry.  At  K  =  10,  a  transverse  vortex  street  appeared,  in  agreement  with  flow 
visualization  experiments. 

The  transverse  vortex  street  observed  at  K  =  12  was  also  reproduced  correctly  and 
for  the  first  time  by  Mostafa  (1987)  using  multi-discrete  vortices.  However,  the  calculated 
forces  were  somewhat  larger  than  those  measured. 

In  later  years,  numerical  experiments  with  co-existing  flows  produced  extremely 
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interesting  flow  features.  For  relative  current  velocities  Vr  (=IVUJ  in  the  range  of  0.7-0.8, 
the  vortices  shed  nearly  symmetrically  at  each  cycle  and  formed  a  three-row  vortex  street  in 
the  range  of  Vr  =  0.6-0.7.  For  Vr  larger  than  about  one,  the  vortex  wake  returned  to  the 
asymmetric  mode,  as  in  a  regular  vortex  street.  Although  the  scope  of  this  work  was  limited 
to  relatively  low  Reynolds  numbers,  the  calculations  of  resistance  in  co-existing  flows 
showed  that  both  the  inertia  and  the  drag  coefficients  for  K  =  4-6  were  in  reasonable 
agreement  with  experimental  data  (Sarpkaya  et  al.  1992). 

Square  cylinders  have  also  been  the  subject  of  numerical  experimentation.  Earlier 
attempts  were  flawed  by  central  differencing  at  large  cell  Reynolds  numbers,  which  led  to 
spatial  oscillations  ahead  of  the  rectangle.  Improvements  were  made  by  using  time 
differencing  and  third-order  upwinding  numerical  schemes  on  the  convective  terms.  The 
technique  just  described  was  tried  for  Reynolds  numbers  under  3,000  (Davis  &  Moore  1982). 
In  1993,  Kato  and  Launder  used  modified  K-epsilon  (k-e)  and  K-omega  (k-Q)  turbulence 
models  on  square  cylinders,  further  improving  the  numerical  results  and  obtaining  a 
markedly  superior  behavior  in  the  near-field  region. 

Aerodynamic  research  on  airfoils  further  demonstrated  the  potential  of  the  k-e  model. 
A  simulation  developed  by  Rogers  (1994),  although  requiring  more  grid  points  than  previous 
ones  (such  as  the  Baldwin-Barth  model),  was  significantly  better  at  computing  maximum  lift 
conditions  and  flap  boundary-layer  separation. 

Whereas  most  of  the  early  numerical  codes  were  generated  by  higher  learning 
institutions  or  government  research  agencies  (such  as  NASA),  a  number  of  very  versatile 
computational  fluid  dynamics  (CFD)  software,  suitable  for  design  and  analysis  purposes, 
is  currently  being  produced  by  the  private  sector.  The  code  utilized  for  this  investigation  is 
a  general-purpose  CFD  code  with  multi-domain  solution  capability  issued  by  CFDRC. 

Initial  results  in  applications  germane  to  the  subject  of  the  current  investigation 
indicated  that  the  program  was  capable  of  predicting  Strouhal  numbers  and  forces  for 
uniform  flows  at  higher  Reynolds  numbers  but  questioned  its  ability  to  capture  the  high 
turbulence  intensity  levels  present  in  the  near- wake  region  (Singhal  &  Awa  1994).  The 
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program  was  further  tested  in  uniform  flow  using  the  k-e  and  renormalization  group  (RNG) 
models  at  Reynolds  numbers  over  106.  Although  it  was  able  to  predict  Strouhal  numbers 
reasonably  well,  the  program  was  not  tested  for  force  coefficients  under  these  flow 
conditions  (Habchi  &  Hufford  1995). 

In  the  current  investigation,  the  CFDRC  software  has  been  used  to  predict  the  forces 
acting  on  circular  cylinders  immersed  in  various  types  of  time-dependent  turbulent  flows. 
It  is  felt  that  the  data  presented  herein  provide  a  good  indication  of  the  code's  performance 
in  this  application. 
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III.  NUMERICAL  REPRESENTATION 


A.  COMPUTATIONAL  METHOD 

The  CFD-ACE  code  simulates  fluid  flow  by  solving  the  partial  differential  equations 
(PDE)  which  govern  the  transport  of  flow  quantities.  Since  the  CFD-ACE  theory  manual 
discusses  the  solution  methodology  at  length,  only  a  brief  outline  of  this  topic  will  be 
presented  here. 

The  governing  equations  for  the  flow  conditions  investigated  are  the  continuity  and 
Navier-Stokes  (NS)  equations.  For  turbulent  flows,  the  code  applies  Favre  (density) 
averaging  to  the  NS  equations,  so  that  each  flow  quantity  4>  is  decomposed  into  a  mean  and 
a  fluctuating  part  according  to: 


4*44"; 


where: 


[1] 


Note  that  overbar  denotes  Reynolds  (time)  averaging,  while  tilde  denotes  Favre  averaging. 
Applying  the  Favre  averaging  procedure  to  the  continuity  and  NS  equations  yields: 


[2] 


and: 


at  1  a  J  1  dx.  dx.  dx.  dx.  3  dx  9  dx  1 


[3] 
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respectively.  By  applying  the  generalized  Boussinesq  eddy  viscosity  concept,  the  Reynolds 
stresses  can  be  treated  as  a  linear  function  of  the  mean  strain  rate,  so  the  Favre  averaged  NS 
(FANS)  equation  can  be  expressed  as: 


J  1  J  J  J  m  i 


[4] 


in  which  r)t  represents  the  turbulent  eddy  viscosity  and  k  is  half  the  trace  of  the  Reynolds 
stress  tensor. 

In  general,  except  for  the  continuity,  each  governing  equation  can  be  expressed  in  a 
common  form  comprised  of  a  transient  term,  a  convective  term,  a  diffusive  term  and  a 
source,  not  all  of  which  will  be  present  at  all  times. 

The  equations  corresponding  to  the  flow  quantities  being  analyzed  are  discretized  by 
the  code  over  finite  control  volumes  and  numerically  integrated,  thereby  yielding  a  set  of 
finite  difference  equations  (FDE).  Although  algebraic,  these  FDEs  are  not  linear  and  are 
therefore  solved  using  an  iterative  process. 

As  mentioned  above,  there  is  no  governing  PDE  available  for  pressure.  The  task  of 
solving  for  this  flow  quantity  while  satisfying  the  continuity  equation  is  accomplished  by  a 
variation  of  the  "Semi-Implicit  Method  for  Pressure-Linked  Equations,  Consistent"  or 
SIMPLEC  algorithm. 

Briefly,  the  SIMPLEC  procedure  can  be  summarized  as  follows: 

1.  Estimate  the  pressure  field. 

2.  Solve  for  the  remaining  flow  quantities  and  check  whether  the  continuity  equation 
is  satisfied.  Since  this  will  not  generally  be  the  case,  an  improved  solution  will  be 
required. 

3.  Find  correction  factors  for  velocities  and  density.  These  can  be  used  to  find  a 
correction  factor  for  pressure. 
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4.  Using  the  correction  factor  obtained  in  the  preceding  step,  re-evaluate  the  pressure 

and  velocity  fields. 

5.  Solve  the  discretized  equations  for  other  flow  variables,  including  turbulence 

quantities,  as  required  by  the  specific  problem. 

6.  Repeat  steps  2  through  5  until  the  solution  converges. 

The  aforementioned  variation  used  in  the  CFD-ACE  code  consists  of  repeating  steps 
3  and  4  a  number  of  times,  thus  updating  the  pressure,  velocity  and  density  fields  prior  to 
proceeding.  These  intermediate  continuity  or  pressure-correction  iterations  have  been  found 
to  enhance  overall  convergence  for  most  flow  problems. 

B.  GRID  GENERATION 

The  grid  used  to  analyze  most  models  is  shown  in  Figure  1.  It  was  developed  using 
the  CFD-GEOM  program  and  is  a  1  meter  by  1  meter  square  subdivided  into  four  domains. 
A  circle  of  diameter  0.05  meters,  located  at  the  geometric  center  of  the  grid,  represents  a 
cylinder  of  infinite  length  about  which  the  fluid  flows.  The  internal  domain  boundaries  are 
located  along  the  diagonals  of  the  square.  Sixty-five  lines  project  radially  from  each  of  the 
four  quadrants  of  the  circle  and  meet  the  sides  of  the  square  at  equispaced  mesh  points. 
Intersecting  the  radial  lines,  one-hundred  and  thirty  other  lines  join  mesh  points  located 
along  the  diagonals  and  spaced  such  that  the  grid  is  much  finer  in  the  vicinity  of  the  circle 
than  in  the  outer  regions. 

Several  requirements  were  considered  during  the  grid  selection  process.  First,  the 
grid  must  be  sufficiently  extensive  to  prevent  any  vortices,  still  strong  enough  to  influence 
the  flow,  from  migrating  beyond  the  boundaries.  Also,  while  it  was  recognized  that 
increasing  the  number  of  cells  (as  well  as  shortening  the  time  interval  between  iterations) 
would  result  in  a  more  accurate  vorticity  computation,  efficiency  considerations  dictated  a 
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finer  mesh  only  in  the  vicinity  of  the  cylinder  walls,  where  the  effect  of  vorticity  is 
strongest.  Unfortunately,  the  vortices  that  separate  and  are  swept  away  diffuse  rapidly  as 
the  mesh  coarsens,  resulting  in  returning  vortices  which  are  probably  weaker  in  the  model 
than  they  would  be  in  actuality. 

When  the  number  of  grid  cells  was  quadrupled,  marginal  improvements  in  the  results 
were  obtained  at  the  cost  of  a  four-fold  increase  in  CPU  time.  It  was  felt,  therefore,  that  the 
grid  just  described  represented  the  best  compromise  between  accuracy  and  computation 
efficiency. 

C.  BOUNDARY  CONDITIONS 

In  the  CFD-ACE  program,  boundary  conditions  can  be  specified  via  drop-down 
menus  or  by  editing  the  file  in  which  all  flow  parameters  are  stored.  Since  this  file  has  an 
".in"  suffix,  it  will  hereinafter  be  referred  to  as  in-file. 

Four  types  of  boundary  conditions  are  used  in  all  the  simulations  performed.  The  top 
and  bottom  sides  of  the  flow  region  are  designated  as  symmetry  boundaries.  This  condition 
implies  a  zero  normal  velocity  and  a  vanishing  gradient  normal  to  the  boundary  for  all 
variables. 

Each  cylinder  quadrant  is  designated  as  a  wall,  implying,  again,  a  zero  normal 
velocity  along  the  boundary. 

The  vertical  sides  of  the  square,  referred  to  as  east  and  west  boundaries,  are  both 
designated  as  exits  with  prescribed  pressure,  resulting  in  a  paradoxical  flow  region  with  two 
outlets  and  no  inlet.  This  apparent  anomaly,  present  in  all  time-dependent  flow  simulations 
performed,  can  be  explained  by  pointing  out  a  peculiarity  of  the  code  nomenclature:  while 
boundaries  designated  as  inlets  do  not  allow  bidirectional  flow,  exits  do.  Several  other  inlet 
and  exit  conditions  may  be  specified  to  accomodate  supersonic,  compressible  and  other  types 
of  flow. 

It  should  be  noted  that  the  flow  conditions  need  not  be  the  same  along  the  entire 
length  of  each  boundary.  It  is  possible,  for  instance,  to  specify  a  higher  flow  velocity  at 
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the  northern  half  of  an  inlet  or  exit  and  a  lower  one  along  the  southern  half.  This  technique 
was  found  particularly  useful  when  flow  perturbations  were  desired  to  precipitate  the 
asymmetrical  shedding  of  vortices  in  uniform  flow. 

A  fourth  type  of  boundary,  called  interface,  is  utilized  to  designate  lines  along  which 
one  domain  meets  with  another  within  the  flow  region. 

D.  FLOW  CONTROL 

In  the  CFD-ACE  program,  velocity  at  an  inlet  or  exit  boundary  can  be  specified 
either  as  a  constant  or  as  one  of  a  limited  number  of  time-dependent  functions.  Pull-down 
menus  can  be  used  to  specify  constant  flows  while  time-dependent  flows  must  be  specified 
by  editing  the  in-file. 

In  the  simulations  performed,  it  was  desired  for  the  ambient  velocity  to  vary  as: 


m  -  u0 


U  cos— t 


[5] 


in  which  U0  represents  a  steady  component,  i.e.,  a  current,  and  Um  the  magnitude  of  the 
fluctuating  component.  Unfortunately,  velocity  was  noted  not  to  follow  harmonic  functions 
whenever  such  a  function  was  invoked.  CFDRC  personnel  acknowledged  this  anomaly 
although  its  underlying  causes  were  not  discussed. 

Flow  can  also  be  controlled  by  making  use  of  the  following  relationship  between 
pressure  and  velocity: 


'I  ■  Tm)  [6] 


and  by  prescribing  a  time-dependent  pressure  at  the  exits.  It  can  be  shown  that,  if  velocity 
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is  to  vary  cosinusoidally,  the  amplitude  of  the  differential  pressure  should  be: 

APm  -  ypU”Ax  [7] 

where  Ax  is  half  the  length  of  the  flow  region.  The  amplitudes  specified  at  the  east  and  west 
exits  should  be  equal  in  magnitude  and  opposite  in  sign.  Since  pressure  and  velocity  are 
off  set  by  id2  radians  by  virtue  of  Equation  [6],  differential  pressure  should  vaiy  sinusoidally, 
resulting  m  a  pressure  which  varies  linearly  in  space  while  its  slope  varies  harmonically  in 

time. 

The  technique  just  described  allows  the  user  to  calculate  the  magnitude  of  the 
oscillating  differential  pressure  required  to  generate  the  desired  flow  velocity.  Whenever  a 
current  is  present,  however,  a  steady  differential  pressure  is  also  required.  It  was  empirically 
determined  that  this  steady  component  is  defined  as: 


APo  “  \  p  Urn  Vr  C  >  where  Vr  =  £8] 

in  which  C  is  dependent  on  K  and  on  the  time  step  used.  The  manner  in  which  this  constant 
was  evaluated,  as  well  as  some  recommended  values  for  different  K  ranges  and  a  time  step 
of  0.02  seconds  are  discussed  in  Appendix  B. 

E.  DATA  REDUCTION  AND  CALCULATION  OF  FORCE  COEFFICIENTS 

Once  the  program  algorithm  has  converged  on  a  set  of  numerical  values  for  the  flow 
quantities  sought  at  each  time  step,  the  results  are  recorded  in  a  file  whose  suffix  is  ".out" 
and  will  therefore  be  referred  to  as  the  out-file. 

In  order  to  reduce  the  data  to  a  useful  format,  it  was  found  necessary  to  transfer  the 
numerical  results  pertinent  to  the  investigation  from  the  out-file  to  a  new  file.  Since  the 
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code  calculates  shear  and  pressure  forces  separately  at  each  boundary  designated  as  a  wall, 
it  is  helpful  if  in  the  course  of  data  transfer  the  contributions  from  each  wall  were 
consolidated  into  a  single  value  for  each  cartesian  axis. 

Once  the  total  horizontal  and  vertical  forces  are  available,  in-line  force  and  lift 
coefficients  are  obtained  as  follows: 


C 
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and: 
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It  is  also  possible  for  the  code  to  record  a  given  flow  quantity  at  a  specified  point  and 
tabulate  the  results  in  a  monitor  file.  Both  flow  quantity  and  location  to  be  monitored  must 
be  specified  by  the  user.  This  capability  was  used  extensively  to  record  velocities  in  order 
to  verify  and  later  fine-tune  the  flow  control  technique  described  in  an  earlier  section. 

Once  the  desired  coefficient  or  flow  quantity  is  tabulated,  its  trend  can  be  plotted  and 
further  analyzed  using  one  of  several  commercially  available  codes. 


F.  INPUT  OPTIONS 

The  CFD-ACE  Command  Language  Manual  describes  the  basic  techniques  used  to 
specify  the  code  inputs  and  the  other  computational  options  available.  Although  drop-down 
menus  can  be  used  in  most  cases,  certain  types  of  inputs  can  only  be  specified  by  editing  the 
in-file.  This  section  will  concentrate  on  specific  aspects  of  input  definition  which  were 
found  particularly  relevant  and  in  some  cases  critical  to  the  type  of  flow  conditions 


13 


investigated.  A  sample  in-file  with  useful  comment  is  presented  in  Appendix  C. 

A  very  important  aspect  of  flow  control  is  the  proper  definition  of  the  coefficients  in 
the  harmonic  functions  which  determine  the  pressure  differential  across  the  flow  region.  Of 
these  coefficients,  Cl  represents  the  steady  state  component  (Ap0in  Equation  [8]),  C2  the 
amplitude  of  the  sinusoidally  oscillating  component  (  Apm  in  Equation  [7]),  and  C3  its 
frequency  in  radians  per  second.  The  remaining  coefficients,  C4  and  C5,  were  set  to  zero  as 
they  are  used  for  cosinusoidal  components,  not  present  in  any  of  the  cases  examined. 

The  user  must  specify  an  initial  and  a  final  time  for  all  time-dependent  flows.  The 
time  increment  is  indirectly  specified  by  inputting  the  number  of  time  steps  to  be  used  and 
should  be  chosen  taking  into  account  the  fineness  of  the  grid  mesh  and  the  flow  velocity. 
In  general  it  is  not  desireable  to  have  the  fluid  travel  a  distance  greater  than  half  of  a  cell 
during  a  single  time  step.  Luckily,  a  modicum  of  experimentation  can  quickly  reveal 
whether  the  improvement  in  output  quality  is  worth  the  longer  computation  time  required 
when  the  time  step  is  shortened. 

It  should  be  noted  that  it  may  be  desirable  for  the  initial  time  to  differ  from  zero 
because  the  program  has  the  capability  of  restarting  from  a  set  of  initial  conditions  recorded 
at  the  end  of  a  previous  run,  allowing  for  two  sets  of  output  data  to  be  spliced  together  and 
plotted  continuously  versus  time. 

The  solution  control  part  of  the  in-file  is  one  where  the  available  options  must  be 
chosen  judiciously  and  some  trial  and  error  may  be  required.  In  most  cases,  the  default 
upwind  spacial  differencing  scheme  and  Euler  temporal  scheme  will  work  satisfactorily. 
Some  preliminary  testing  of  other  options  led  to  the  observations  which  follow. 

Use  of  the  Crank-Nicholson  scheme  cannot  be  prescribed  using  the  appropriate  drop¬ 
down  menu  as  this  will  result  in  a  syntax  error  within  the  in-file.  In  general,  however,  this 
technique  does  not  produce  appreciable  changes  with  respect  to  the  results  achieved  using 
the  default  Euler  mehod. 

While  the  upwind  spacial  differencing  method  is  very  robust,  central  differencing 
may  become  unstable,  especially  when  used  in  conjunction  with  fine  grid  meshes.  This 
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problem  can  usually  be  overcome  by  the  appropriate  selection  of  a  blending  factor:  in  most 
models  in  which  central  differencing  was  used,  a  blending  factor  of  0.3  was  specified;  0.5 
was  successfully  substituted  in  models  which  diverged.  The  greater  the  blending  factor,  the 
more  the  techinque  will  have  an  upwind  character. 

Other  spacial  differencing  schemes  did  not  work  satisfactorily  on  the  type  of 
computer  used.  A  later  version  of  the  program  corrected  this. 

The  number  of  solution  iterations  and  sub-iterations  for  each  flow  quantity  can  only 
be  selected  using  some  experimentation  or  previous  experience  and  it  will  differ  for  each 
type  of  problem  analyzed.  In  general,  all  iterative  processes  should  be  shortened  as  much 
as  possible  without  compromising  accuracy.  This  can  only  be  done  by  solving  a  problem 
with  a  known  solution  and  changing  the  number  of  iterations  until  a  satisfactory  result  is 
obtained.  The  program  allows  for  visualization  of  the  trend  in  the  residuals  of  each  flow 
quantity  calculation,  so  the  slope  of  residual  graphs  is  also  helpul  in  determining  the 
optimum  number  of  iterations. 

The  code  offers  a  choice  of  six  different  turbulence  models,  only  three  of  which  were 
utilized  in  the  course  of  this  investigation.  The  standard  k-e  model  was  found  to  work 
satisfactorily  in  most  circumstances.  Very  minor  changes  were  noted  in  the  output  when  the 
RNG  model  was  used.  Whenever  the  grid-flow  velocity  combination  is  such  that  y+  (=(y/v 
)(t  /p  )0  5)  is  smaller  than  about  1 1  at  the  grid  cell  center  closest  to  the  cylinder  wall,  use  of 
the  Low  Reynolds  model  is  appropriate  and  results  in  better  numerical  accuracy.  Since  in 
time-dependent  flows  velocity  varies,  the  user  needs  to  determine  whether  the  values  of  y+ 
warrants  selecting  the  Low  Reynolds  model.  All  the  above  models  require  the  user  to  input 
two  parameters,  turbulence  kinetic  energy  and  turbulence  length  scale.  It  was  found  that 
varying  the  former  from  a  value  corresponding  to  3%  to  10%  turbulence  caused  minimal 
changes  in  the  results,  so  3%  was  used  throughout.  The  latter  was  varied  between  0.09  D/2 
and  0.09  D  (D  being  the  cylinder  diameter)  without  noting  significant  differences. 
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IV.  DISCUSSION  OF  RESULTS 


A.  INTRODUCTION 

The  results  will  be  presented  in  the  following  order.  First,  a  discussion  of  the  grid 
and  other  code  inputs  used  in  the  calculations.  Second,  a  detailed  discussion  of  the  in-line 
and  transverse  force  coefficients  as  functions  of  the  selected  governing  parameters.  Finally, 
a  comparison,  whenever  possible,  of  the  calculated  and  experimental  values.  This  section 
is  then  completed  with  plots  of  vorticity  and  streamlines  for  representative  values  of  flow 
parameters.  Throughout,  a  special  attempt  is  made  to  concentrate  on  the  most  informative 
range  of  the  calculations  rather  than  on  a  massive  presentation  of  the  data. 

B.  GRID  AND  INPUT  PARAMETERS 

As  discussed  in  Section  III  B,  the  code  used  allowed  one  to  generate  various  types 
of  grids.  For  the  problem  under  consideration,  the  axial  symmetry  of  the  cylinder,  the 
problem's  two-dimensionality  and  oscillating  nature  (left  and  right  boundaries),  and  the 
extreme  importance  of  the  accurate  calculation  of  the  velocity  gradients  at  and  near  the 
cylinder  surface  were  kept  in  mind.  Several  trials  and  tribulations  finally  resulted  in  the  grid 
shown  in  Fig.  1 .  It  consists  of  four  domains.  The  smallest  cell  size  is  3.05*  1 0‘5  D  along  the 
circumference  of  the  cylinder  and  2.56*1  O'5  D  in  the  radial  direction.  The  largest  cell,  at  the 
outer  boundaries,  is  0.313  D  by  0.245  D.  The  size  of  the  computational  domain  is  20  D  by 
20  D. 

Careful  consideration  had  to  be  given  to  the  selection  of  the  grid  size  for  reasons 
beyond  the  accurate  calculation  of  the  velocity  gradient  on  the  cylinder.  The  fact  that  the 
flow  oscillates  makes  the  vortices  generated  during  one  half-cycle  return  toward  their 
inception  point.  Thus,  if  a  vortex  does  not  retain  its  integrity  as  it  transits  areas  of  coarser 
mesh,  then  it  will  return  to  the  cylinder  with  artificial  diffusion  above  and  beyond  that 
imposed  by  the  prevailing  turbulence.  Furthermore,  the  stability  of  the  flow,  regardless  of 
where  the  cell  may  be,  is  dictated,  among  other  parameters,  by  the  Von  Neumann  criterion. 
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It  states  that  flow  at  a  given  point  should  not  displace  more  than  one  half  of  the  cell  on  which 
the  fluid  particle  resides.  This  is  a  difficult  condition  to  satisfy  even  in  unidirectional  steady 
flows.  Since  in  periodic  flows,  with  or  without  current,  the  magnitude  of  the  ambient 
velocity  fluctuates  between  zero  and  a  prescribed  maximum,  then  one  either  allows  violation 
of  the  Von  Neumann  criterion  at  times  of  peak  velocity,  so  that  global  instability  is  not 
allowed  to  dominate  the  flow,  or  is  penalized  in  terms  of  grid  fineness,  and  consequently 
CPU  time,  by  satisfying  the  criterion  throughout.  The  experience  of  others,  as  well  as  the 
present  one  have  demonstrated  that  the  latter  is  impractical. 

In  the  present  simulations,  a  compromise  was  made  to  maintain  sufficient  accuracy 
while  sacrificing  neither  physics  nor  a  reasonable  CPU  time.  It  is  as  a  result  of  the  foregoing 
considerations  that  the  quantity  UmAt/£,  where  £  is  the  grid  cell's  characteristic  length,  was 
allowed  to  reach  values  on  the  order  of  50  in  the  vicinity  of  the  cylinder  walls.  It  should  be 
noted,  of  course,  that  the  quantity  u(t)At/£,  i.e.  the  ratio  of  the  actual  distance  traveled  by  the 
fluid  to  the  computational  cell's  characteristic  length,  is  expected  to  be  substantially  less  than 
UmAt/£  in  the  immediate  vicinity  of  the  solid  wall  boundary. 

In  the  results  to  be  presented,  the  flow  parameters  chosen  were  the  Keulegan- 
Carpenter  number  K,  the  relative  velocity  parameter  Vr  and  the  frequency  parameter  p  or 
the  Reynolds  number  Re.  The  majority  of  the  calculations  were  performed  with  a  time 
interval  At  of  0.02  seconds.  As  will  be  noted  in  connection  with  the  sensitivity  analysis 
discussion,  the  time  interval,  and  other  code  and  problem  related  parameters,  such  as  the 
turbulence  models,  were  judiciously  varied.  It  is  fully  realized  that  a  numerical  simulation, 
like  a  physical  experiment,  produces  only  a  data  point  on  the  basis  of  the  input  parameters 
chosen  and  discretization  schemes  employed  in  the  code.  Thus,  the  insight  that  one  develops 
through  familiarity  with  the  phenomenon  is  brought  to  bear  on  the  selection  of  the  input 
parameter  and  on  the  comparison  of  the  results  with  those  obtained  experimentally.  While 
the  calculations  undoubtedly  suffered  from  boundary  constraints  and  artificial  diffusion,  the 
experiments  have  their  own  corresponding  limitations  with  more  or  less  uncertainties. 

Once  the  parameters  were  chosen,  the  code  was  capable  of  producing  the  in-line  and 
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transverse  forces  and  the  necessary  information  for  the  plot  of  vorticity  and  streamlines.  The 
calculations  have  been  carried  out  for  a  minimum  of  seven  and  a  maximum  of  30  cycles. 

C.  IN-LINE  AND  TRANSVERSE  FORCE  COEFFICIENTS 

These  will  be  discussed  for  p  =  1772,  2487,  and  4122.  In  case  of  p  =  1772,  the 
analysis  was  carried  out  for  six  Vr  values,  corresponding  to  six  K  values  ranging  from  8.6 
to  36.  For  p  =  2487,  for  five  Vr  values  corresponding  to  five  K  values  ranging  from  10.5  to 
46.  Finally,  for  p  =  4122,  for  three  values  of  Vp  corresponding  to  three  values  of  K  ranging 
from  5.2  to  24.  In  addition,  for  p  =  1772  only,  the  parameter  Vk  (=  KVr  =  U0T/D)  was 
assigned  three  different  values  in  order  to  compare  the  results  with  those  obtained 
experimentally  wherever  possible.  For  the  other  two  p  values,  only  single  values  of  Vk  were 
used.  In  fact,  it  is  because  of  this  reason  that  there  was  a  one-to-one  correspondence  between 
Vr  and  K  for  each  combination  of  p  andVk. 

Figures  2a  through  2d  show  the  in-line  force  coefficient  (Cn)  and  the  transverse  force 
coefficient  (Ctl)  for  p=1772,  K=8.6  and  Vk=2.05.  The  basic  difference  between  the  first  two 
and  the  second  two  figures  is  that  the  former  were  produced  using  the  k-e  turbulence  model 
and  the  latter  using  the  low  Reynolds  number  (Low-Re)  model.  The  purpose  of  this  exercise 
was  to  demonstrate  the  consequences  of  the  use  of  two  turbulence  models. 

It  is  realized  that,  for  the  case  being  discussed,  the  maximum  Re  is  approximately 
15,000  and  that,  during  any  given  cycle,  the  flow  undergoes  a  transition  from  laminar  to 
turbulent.  Thus,  the  use  of  a  turbulence  model  such  as  the  k-e,  tailored  for  use  in  high  Re 
steady  flows,  may  not  be  entirely  appropriate  for  an  unsteady  flow  at  relatively  low  Re.  A 
comparison  of  Figs.  2a  and  2c  shows  that,  whereas  in  the  former  Cn  does  not  have  an 
inflection  point  even  as  late  as  the  fourth  cycle,  in  the  latter  an  inflection  point  is  evident 
shortly  after  the  first  cycle.  The  fact  that  the  evolution  of  the  inflection  point  is  directly 
related  to  the  occurrence  of  asymmetric  vortex  shedding  is  an  indication  of  the  evolution  of 
the  lift  force.  In  fact,  a  comparison  of  Figs.2b  and  2d  shows  that  there  is  no  significant  lift 
in  the  former  relative  to  the  latter.  In  other  words,  the  use  of  the  Low-Re  model,  thought  to 
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be  more  appropriate  for  low  Re  flows,  gives  rise  to  transverse  force  and  vortex  shedding  after 
a  few  cycles  from  the  inception  of  flow. 

A  comparison  of  the  turbulence  models  cannot  be  made  for  all  the  cases  considered 
in  this  study  because  of  obvious  time  limitations.  Nevertheless,  an  effort  is  made  to  combine 
such  sensitivity  analyses  into  a  section  to  be  presented  later. 

Returning  to  the  presentation  of  the  results,  and  in  particular  to  the  case  of  p  =  1772, 
Figs.  3  and  4  show  the  Q,  for  Vk=2.05  for  K=15  and  36.5.  It  is  evident  from  a  brief  perusal 
of  these  figures  that  the  in-line  force  assumes  a  quasi-steady  state  within  four  to  six  cycles. 
The  shape  of  the  C„  graph,  particularly  in  Fig.  4,  is  reminiscent  of  an  odd  harmonic  function, 
with  the  inflection  points  near  the  upper  left  and  the  lower  right  of  each  cycle  associated,  as 
noted  earlier,  with  the  alternate  shedding  of  the  vortices  from  the  top  and  bottom  of  the 
cylinder. 

Figures  5a,  6a,  6c  and  7  show  CH  again  for  p=1772  for  Vk=4  and  for  three  different 
K  values.  The  most  noteworthy  among  the  four  plots  is  Fig.  6a,  where  Cu  exhibits  irregular 
behavior  after  approximately  eight  cycles,  having  earlier  reached  a  quasi-steady  state.  This 
seemingly  strange  occurrence  is  due  to  the  phenomenon  known  as  the  transverse  vortex 
street.  For  K  values  in  the  range  of  10<k<13,  the  vortices  shed  from  the  cylinder  move  in 
the  transverse  direction  and  only  on  one  side  of  the  cylinder.  This  means  that  the  vortices 
shed  in  each  cycle  go  around  the  cylinder  and  find  their  way  in  the  positive  or  negative 
vertical  direction,  thereby  creating  asymmetrical  and  unsteady  in-line  forces.  This 
phenomenon  has  been  discussed  in  detail  by  Sarpkaya  and  Isaacson  (Sarpkaya  and  Isaacson, 
1981).  In  fact.  Fig.  6b  shows  that  a  transverse  force  develops  after  approximately  seven 
cycles  and  reflects  the  shedding  of  the  vortices  at  the  top  and  bottom  of  the  cylinder  for  a 
number  of  cycles.  It  should  be  noted  that  Figs.  6a  and  6b  were  developed  using  the 
renormalization  group  theory  (RNG)  turbulence  model,  whereas  the  low-Re  model  was  used 
to  generate  Figs.  6c  and  6d.  Even  though  the  Cn's  do  not  differ  significantly,  the  transverse 
force  coefficients  (C,,)  show  the  dramatic  effect  of  the  choice  of  turbulence  model  used  in  the 
evolution  of  the  vortex  asymmetry.  As  will  be  noted  later,  not  all  predictions  of  different 
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turbulence  models  will  yield  the  same  result  even  when  the  same  set  of  governing  parameters 
is  used. 

Figures  8  and  9a  show  Q,,  again  for  p=l  772,  for  a  higher  value  of  Vk  and  for  K=1 0.3 
and  36.5.  Figure  9a  is  particularly  interesting  in  that  Cn  is  asymmetrical  in  shape  and 
magnitude  even  though  K  is  as  large  as  36.5.  Normally,  one  would  expect  a  steady  flow  like 
behavior  for  Vk=0  and  relatively  large  K  values.  An  examination  of  Figs.  9a  and  9b  indicates 
that  neither  Cfl  nor  Ctl  are  symmetrical  in  the  presence  of  a  current,  whereas  Figs.  9c  and  9d, 
which  correspond  to  flows  with  no  current  but  otherwise  at  identical  values  of  K  and  p,  show 
that  both  Cn  and  Ctl  are  as  symmetrical  as  they  can  be  after  an  initial  stage  of  flow 
development.  Thus,  the  presence  of  current  has  a  significant  impact  on  the  topology  of 
vortex  shedding  although  the  superposition  of  collinear  waves  and  current  do  not  necessarily 
lead  to  periodically  alternating  vortex  shedding.  It  is  because  of  this  reason  that  Cn  in  Fig. 
9a  exhibits  symmetry  neither  in  magnitude  nor  in  shape. 

The  next  case  to  be  considered  is  for  an  intermediate  p  value  of  2487.  Figures  10 
through  12  show  Cu  for  Vk=5.01  and  for  three  values  of  K.  The  most  significant  aspect  of 
the  graphs  is  that  the  amplitude  of  the  oscillation  decreases  with  increasing  K  and  the  shape 
of  the  Cn  trace  becomes  less  sinusoidal,  reflecting,  at  each  half  cycle,  the  asymmetric 
shedding  of  the  vortices. 

Figure  13  shows  the  Ca  for  the  highest  p  value  dealt  with  in  this  investigation  for 
Vk=2.07.  The  force  oscillations  are,  as  expected,  fairly  sinusoidal  since,  for  small  K  values, 
one  does  not  expect  significant  vortex  shedding. 

Figures  14a  and  14b  for  Vk=2.07  and  (3=4122  shows  what  happens  to  in-line  and  lift 
forces  in  the  most  dramatic  range  of  K.  When  K  is  in  the  range  between  10  and  20,  vortex 
shedding  becomes  very  chaotic,  Cu  decreases  significantly  and  Cn  drops  down  relative  to  the 
case  of  K=5. 

The  presence  of  current  further  complicates  the  flow  topology  with  respect  to  the  no¬ 
current  case  in  this  most  sensitive  region  of  the  K  values.  The  flow  does  not  have  any  sort 
of  symmetry  about  either  axis.  Even  though  the  prediction  of  forces  in  this  region  is  rather 
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difficult,  there  is  one  beneficial  effect  deriving  from  the  smallness  and  the  chaotic 
character  of  the  lift,  i.e.  the  minimization  or  practical  elimination  of  the  tendency  of  the 
cylinder  to  develop  a  dynamic  response  or  hydroelastic  oscillations.  However,  for  K=24 
(P=4122  and  Vk=2.07),  the  vortex  shedding  becomes  quasi-periodic  and  the  C„  increases 
considerably  as  shown  in  Figs.  15a  and  15b. 

A  harmonic  analysis  of  the  lift  coefficient  has  shown  that  the  ratio  of  the 
frequency  of  vortex  shedding  to  the  flow  oscillation  frequency  is  in  the  range  of  3  to  4 
depending  on  Vk.  Such  an  analysis  is  not  shown  here  but  may  be  found  in  Storm's  work. 
(Storm,  1984). 

D.  FLOW  VISUALIZATION 

Figures  16a  through  16f  show  the  evolution  of  vorticity  for  K=8.78  and  P=1772. 
In  this  particular  flow  regime,  vortices  are  formed  and  shed  from  the  body  in  quasi- 
symmetrical  pairs,  so  that  no  significant  amount  of  lift  is  generated.  This  observation  is 
reinforced  by  examining  Figs.  17a  through  17f,  which  show  the  streamlines  corresponding 
to  the  same  flow  regime. 

Figures  18a  through  18f  show  the  more  complex  flow  topology  resulting  from  a 
slightly  higher  K  and  a  much  higher  p.  (10.5  and  2487  respectively).  The  corresponding 
streamlines,  shown  in  Figs.  19a  through  19f,  clearly  indicate  that  the  vortex  pattern  is 
asymmetrical,  so  that  a  sizable  lift  force  results. 

E.  NUMERICAL  RESULTS 

The  values  of  the  maximum  in-line  force  coefficients  and  the  corresponding 
governing  parameters  are  listed  in  Tables  1  through  5.  Unfortunately,  experimental  results 
were  available  only  for  one  set  of  the  governing  parameters.  The  general  trends  noted  are 
that  the  agreement  with  both  the  Morison's  equation  and  with  experimental  data  is  better 
for  K  values  less  than  about  10  and  larger  then  20.  This  confirms  that  in  a  K  range  close 
to  that  in  which  a  transverse  vortex  street  can  be  expected  to  develop,  the  flow  topology 
is  particularly  challenging  to  capture.  Another  trend  noted  was  that,  as  the  value  of  Vk 
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increases,  the  agreement  between  the  computed  CH  values  and  those  produced  by  the 
Morison's  equation  will  deteriorate  slightly.  The  availability  of  additional  experimental  data 
at  higher  Vk  would  be  required  to  draw  more  insight  from  the  observed  trend. 

F.  SENSITIVITY  ANALYSIS 

A  fairly  detailed  sensitivity  analysis  for  no-current  cases  was  conducted  by  Hanson 
(Hanson,  1995).  Since  in  that  investigation  the  same  computational  grid,  turbulent 
parameters,  and  time  increment  were  used  as  those  employed  in  this  study,  most  of  the 
findings  are  valid  for  the  cases  examined  herein. 

The  analysis  conducted  within  the  scope  of  this  work  involves  primarily  turbulent 
models  and  time  increments.  Since,  due  to  time  constraints,  this  analysis  was  performed  on 
representative  cases  only,  some  of  the  observations  made  will  be  mostly  qualitative  in  nature. 

In  all  cases  examined,  the  turbulence  parameters  prescribed  corresponded  to  3% 
turbulence  and  an  integral  turbulent  length  scale  of  0.00225.  The  fluid  density  was  998 
Kg/m3  and  the  kinematic  viscosity  1*106  m2/s.  A  spacial  central  differencing  scheme  with 
a  blending  factor  of  0.3  was  used,  along  with  an  Euler  temporal  scheme.  Other 
computational  parameters  can  be  obtained  from  Appendix  C,  which  is  a  fairly  typical  in-file. 

Several  turbulence  models  were  tested  for  the  case  corresponding  to  K=12.19, 
P=1772,  Vk=4.  While  these  all  produced  values  of  Ca  within  3%  of  each  other,  the  k-e 
model  produced  no  inflection  point  after  the  second  oscillation  and  the  vortex  pattern 
appeared  to  be  still  symmetrical  after  7  cycles  of  flow  oscillation.  In  contrast,  the  traces 
produced  by  the  RNG  and  Low  Re  models  are  shown  in  Figs.  6a  through  6d. 

As  previously  stated,  the  quantity  UmAt/<;  was  allowed  to  reach  values  in  excess  of 
50  at  the  cylinder  walls.  The  case  corresponding  to  K=24,  p=4122  (Re=  100,000),  for  which 
the  initial  value  of  UmAt/e  was  62.3  was  examined  as  it  represents  one  of  the  highest  values 
of  that  quantity  encountered  in  this  study.  The  time  increment  was  first  adjusted  so  that 
UmAt/q  was  approximately  12.5.  This  resulted  in  a  reduction  in  C,,  from  1.22  to  1.05.  (A 
similar  procedure  for  K=46  and  p=2487  resulted  in  a  reduction  in  Ca  from  1.1  to  0.8.)  A 
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further  reduction  in  UmAt/?  to  0.623  did  not  result  in  a  sizeable  change  with  respect  to  the 
case  corresponding  to  UmAt/c=12.5  (less  than  1%).  Very  little  effect  on  Cu  (<2%)  was  also 
noted  when  UmAt/?  was  changed  from  approximately  15  to  less  than  0.5  for  a  case  where 
K=12.19  and  p=1772. 

It  is  reasonably  safe,  therefore,  to  conclude  that,  for  the  particular  computational  grid 
used,  values  of  UmAt/g  on  the  order  of  15  should  be  acceptable.  Although  there  is  not 
sufficient  data  to  pinpoint  the  critical  value  of  UmAt/?  for  which  a  significant  impact  on  the 
result  should  be  expected,  this  value  will  almost  certainly  be  between  about  20  and  50. 
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TABLE  1 


Computed  in-line  force  coefficients  for  different  values  of  K  and  p  =  1772, 
Vk  =  U0K  =  2.05 


K 

vr 

Cg 

(calc.) 

Cn  (Morison) 

8.6 

0.22 

2.39 

2.37 

12.1 

0.16 

1.75 

2.12 

15.0 

0.14 

1.49 

1.85 

17.8 

0.12 

1.31 

1.66 

25.8 

0.09 

1.05 

1.30 

36.5 

0.07 

0.92 

1.08 
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TABLE  2 


Computed  and  experimental  in-line  force  coefficients  for  different  values  of  K  and  p=1772 
Vk  =  UoK  =  4.0 


K 

vr 

Q,  (calc.) 

c„ 

(Mor) 

Cn  (Exper) 

8.78 

0.43 

2.39 

2.78 

2.60 

12.19 

0.31 

1.83 

2.42 

2.25 

17.82 

0.23 

1.33 

1.87 

1.70 

26.07 

0.16 

1.09 

1.43 

1.10 

TABLE  3 


Computed  in-line  force  coefficients  for  different  values  of  K  and  p  =  1772, 
Vk  =  U0K  =  6.05 


K 

vr 

Q,  (calc.) 

C;, 

(Mor) 

10.3 

0.56 

2.11 

3.56 

12.1 

0.49 

1.84 

3.17 

25.8 

0.25 

1.14 

1.61 

36.5 

0.18 

1.02 

1.31 
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TABLE  4 


Computed  in-line  force  coefficients  for  different  values  of  K  and  p  =  2487, 
Vk  =  U0K  =5.01 


K 

vr 

Ca  (calc.) 

c„ 

(Mor) 

10.5 

0.47 

2.10 

3.05 

12.4 

0.40 

1.84 

2.56 

17.2 

0.30 

1.48 

1.86 

29 

0.19 

1.18 

1.14 

46 

0.13 

1.10 

0.81 

TABLE  5 


Computed  in-line  force  coefficients  for  different  values  of  K  and  p  =  4122, 
Vk  =  UoK  =  2.07 


K 

vr 

Ca  (calc.) 

Cd 

(Mor) 

5.1 

0.41 

4.14 

3.44 

15 

0.14 

1.58 

1.33 

24 

0.11 

1.22 

0.97 

27 
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V.  CONCLUSIONS 


The  purpose  of  this  investigation  was  a  critical  assessment  of  the  numerical  analyses 
of  two-dimensional,  sinusoidally  oscillating,  turbulent  flows  at  relatively  large  Re,  K  and  Vr 
A  number  of  turbulence  models,  including  the  standard  k-e,  RNG  based  k-e,  and  Low-Re 
models  have  been  employed.  The  emphasis  was  on  the  validation  and  accuracy  of  the  time- 
dependent  computations.  Observations  relating  to  both  numerical  and  physical  validation 
lead  to  the  following  remarks: 


(1)  The  finite  differencing  formulation  of  the  Favre  Averaged  Navier  Stokes 
equation  (FANS)  can  reasonably  solve  flows  within  a  wide  range  of  Re  using  a 
modified  SIMPLEC  method.  However,  the  predicted  forces  will  almost  always 
be  smaller  than  those  obtained  experimentally. 

(2)  The  turbulence  models  used  do  not  fully  capture  vortex  strength  and  prematurely 
dissipate  the  vortices  relative  to  experimental  observations  and  measurements. 
All  models  used  herein  (standard  k-e,  RNG  based  k-e,  and  Low-Re  models)  fail 
to  predict  accurately  the  exact  size  of  the  vortices  and  the  high  turbulence 
intensity  levels  present  in  the  near- wake  of  time-dependent  flows  (subjected  to 
unsteady  pressure  gradients).  This  is  the  primary  cause  of  the  aforementioned 
tendency  to  underpredict  the  force  coefficients. 

(3)  The  use  of  much  finer  grids,  reduced  time  increments  and  higher-order  spacial 
and  temporal  differencing  schemes  (at  the  expense  of  increased  CPU  time)  will 
not  always  increase  the  accuracy  of  the  predictions  since  the  validity  of  the 
turbulence  models  for  flows  subjected  to  extra  strains  such  as  the  time- 
dependent  pressure  gradients  remains  unknown. 

(4)  Turbulence  is  at  present,  and  is  likely  to  remain  for  a  long  time,  the  greatest 
roadblock  to  computational  fluid  dynamics.  Some  appreciation  of  the 
incomplete  knowledge  bases  (both  numerical  and  experimental),  retrofitting  of 
data,  an  the  assessment  of  their  consequences  are  necessary  to  achieve  often 
qualitative  and  occasionally  quantitative  simulations,  particularly  in  time- 
dependent  flows.  This  is  a  compromise  between  expectations  and  achievables 
and  between  physically  relevant  dynamics  and  specific  quantitative  results. 
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Figure  1.  Computational  Grid. 
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In-line  Force  coeff 


(Low  Re  turbulence  model  used) 


Figure  2d.  Transverse  force  coefficient  versus  time.  K=8.6;  P=1772;  Vk=2.05 
(Low  Re  turbulence  model  used) 
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Figure  6b.  Transverse  force  coefficient  versus  time.  K=12.19;  P=1772;  Vk=4. 
Inception  of  asymmetric  vortex  shedding  at  t/T~7.  (RNG  turb.  model) 
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Cycles,  t/T 

Figure  6c.  In-line  force  coefficient  versus  time.  K=12.19;  p=1772;  Vk=4. 
Inception  of  asymmetric  vortex  shedding  at  t/T~3.  (Low  Re  turb.  model) 
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Figure  6d.  Transverse  force  coefficient  versus  time.  K=12.19;  (3=1772;  Vk=4. 
Inception  of  asymmetric  vortex  shedding  at  t/T~3.  (Low  Re  turb.  model) 
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Cycles,  t/T 

Figure  8.  In-line  force  coefficient  versus  time.  K=10.3;  P=1772;  Vk=6.05 
Inception  of  asymmetric  vortex  shedding  at  t/T =4 
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Lift  coeff 


Figure  9a.  In-line  force  coefficient  versus  time.  K=36.5;  p=1772;  Vk=6.05 
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Cycles,  t/T 

Figure  10.  In-line  force  coefficient  versus  time.  K=10.5;  P=2487;  Vk=5.01 
Inception  of  asymmetric  vortex  shedding  at  t/T*5 


Cycles,  t/T 


Figure  13.  In-line  force  coefficient  versus  time.  K=5.1;  p=4122;  Vk=2.07 
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coefficient  versus  time.  K=15;  P=4122;  Vk=2.07 
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APPENDIX  B. 

RECOMMENDED  VALUES  FOR  THE  CURRENT  VELOCITY- 
PRESSURE  DIFFERENTIAL  CORRELATION  COEFFICIENT 

It  was  found  experimentally  that  the  coefficient  C  in  equation  [8]  of  section  HI  D 
assumes  different  values  as  K  and  the  computational  time  interval  change.  Since  a 
considerable  number  of  calculations  were  performed  with  a  time  interval  of  0.02  seconds, 
the  following  table  of  recommended  values  for  C,  valid  for  time  steps  on  the  order  of 
0.02  seconds,  is  included.  It  should  be  realized  that  the  recommended  values  constitute 
only  a  starting  point  and  that  further  fine  tuning  might  be  necessary  depending  on  how 
accurately  the  value  of  U0  is  to  be  duplicated.  One  can  infer  from  the  table  that,  as  K 
increases,  C  becomes  closer  to  0.5.  In  fact,  when  0.5  was  used  in  one  steady  flow  case 
(which  can  be  thought  as  a  flow  with  an  infinite  value  of  K),  the  correct  value  of  ambient 
velocity  was  obtained. 

As  far  as  the  dependence  of  C  on  the  time  interval,  as  a  thumbrule  C  will  assume 
values  closer  to  0.5  as  the  time  interval  is  shortened.  Unfortunately,  the  nature  of  this 
dependence  has  been  noted  only  toward  the  end  of  this  investigation  so  that  not  even  an 
approximate  quantitative  definition  of  it  can  be  formulated  at  this  time. 


Range  of  K  values 

Recommended  C  value 

<10 

0.5892 

10<K<18 

0.5721 

.  18<K<28 

0.5561 

28<K<38 

0.5411 

K>38 

0.5261 
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APPENDIX  C.  SAMPLE  IN-FILE 


TITLE  'K=36.5,  Beta=1772,  0=158.49,  Vk=6.05,  T=1.41  coarse  grid,  no  disturbance* 
* 

* 

* 

* 

* 

**** ** ****** **********  Model  *********************** 

* 

* 

MODEL  tcgll5t 
FUNCTIONS 

PERIODIC  WAVEPTOP  Cl=l  58.49  C2=2874.69  C3=4.453  C4=0  C5=0 
PERIODIC  WAVEPBOT  Cl=l 58.49  C2=2874.69  C3=4.453  C4=0  C5=0 
PERIODIC  WAVENTOP  Cl=-1 58.49  C2=-2874.69  C3=4.453  C4=0  C5=0 
PERIODIC  WAVENBOT  Cl=-1 58.49  C2=-2874.69  C3=4.453  C4=0  C5=0 
END 

* 

GEOMETRY 
GRID  2D  BFC 

READ  GRID  FROM  tcg.PFG 

*  No  Cell  Types  (Blockage  or  Solid)  specified. 

END 

* 

PROBLEM_TYPE 
SOLVE  FLOW  TURBULENCE 
UNSTEADY  TF  =  0  TL  =  10  STEPS  =  500 
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END 


PROPERTIES 

DENSITY  CONSTANT  998 

VISCOSITY  CONSTANT  KINEMATIC  le-06 

END 

* 

MODELS 

TURBULENCE  MODEL  KE 

END 

* 

***  Boundary  Conditions  *** 


BOUND  ARY_CONDITION  S 

*  boundary  condition:  INTERFACE_1_1 
INTERFACE  1  129  1  1  SOUTH 

*  boundary  condition:  wall 
WALL  1  1  1  64  WEST 
U  =  0  V  =  0 

*  boundary  condition:  outritghttop 
EXIT_P  129  129  32  64  EAST 

U  =  0  V  =  0  P  =  FUNCTIONJWAVENTOP  K  =  .000753  D  =0  L  =  0.00225  T 
293 

*  boundary  condition:  outrightbot 
EXITP  129  129  1  31  EAST 

U  =  0  V  =  0  P  =  FUNCTION  WAVENBOT  K  =  .000753  D  =0  L  =  0.00225  T 
293 

*  boundary  condition:  INTERFACE_1_2 
INTERFACE  1  129  64  64  NORTH 
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DOMAIN  2 

*  boundary  condition:  INTERFACE_2_1 
INTERFACE  64  64  1  129  EAST 

*  boundary  condition:  wall 
WALL  1  64  1  1  SOUTH 
U=0 V=0 

*  boundary  condition:  symm 
SYMMETRY  1  64  129  129  NORTH 

*  boundary  condition:  INTERFACE_2_2 
INTERFACE  1  1  1  129  WEST 
DOMAIN  3 

*  boundary  condition:  INTERFACE_3_1 
INTERFACE  1  129  1  1  SOUTH 

*  boundary  condition:  outlefttop 
EXIT_P  1  1  32  64  WEST 

U  =  0V  =  0  P  =  FUN  CTION_W  A  VEPT  OP  K  =  .000753  D  =  0  L  =  0.00225  T 
293 

*  boundary  condition:  outleftbot 
EXIT_P  11131  WEST 

U  =  0V  =  0  P  =  FUNCTION_WAVEPBOT  K  =  .000753  D  =  0  L  =  0.00225  T 
293 

*  boundary  condition:  wall 
WALL  129  129  1  64  EAST 
U  =  0  V  =  0 

*  boundary  condition:  INTERFACE32 
INTERFACE  1  129  64  64  NORTH 
DOMAIN  4 

*  boundary  condition:  INTERFACE_4_1 
INTERFACE  64  64  1  129  EAST 


53 


*  boundary  condition:  symm 
SYMMETRY  1  64  1  1  SOUTH 

*  boundary  condition:  wall 
WALL  1  64  129  129  NORTH 
U  =  0  V  =  0 

*  boundary  condition:  INTERFACE_4_2 
INTERFACE  1  1  1  129  WEST 

*  No  Momentum  Resistances  Defined 

END 

* 

INITXAL_CONDITIONS 

*  Full  field  initial  conditions 

U  =  -1.079  V  =  0P  =  0T  =  293  K  =  .000753  D  =  0  L  =  0.00225 

*  RESTART  FROM  tcg43.2000.AUR 

END 

* 

^^♦^jltjJejlc^jle*************  *********************** 


SOLUTIONCONTROL 
ALGORITHM  SIMPLEC 
S_SCHEME  UPWIND  RHO  K  D 
S_SCHEME  CENTRAL  U  V 
SJBLENDING  0.3  U  V 
T_SCHEME  EULER 
ITERATIONS  10 
CITERATIONS  1 
SOLVER  WHOLE_I  U  V  PP  K  D 
S  ITERATIONS  4  U  V 
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SJTERATIONS  20  PP 
SJTERATIONS  10  K  D 
INERTIALFACTOR  0.2  U  V 
INERTIALF  ACTOR  0.2  K  D 
RELAX  0.4  P 
RELAX  1  RHO  T  VIS 
MINVAL  -le+20  U  V 
MINVAL  -le+20  P 
MINVAL  le-06  RHO 
MINVAL  le-10  T  VIS 
MINVAL  le-10  K  D 
MAXVAL  le+20  U  V 
MAXVAL  le+20  P  RHO 
MAXVAL  5000  T 
MAXVAL  le+20  VIS 
MAXVAL  le+20  K  D 

END 

* 

OUTPUT 

P_FORCE  ON  PMIN  =  -l.e+05 
PLOT3D  ON  FORMATTED 
SCALAR  FELE  1  RHO  P  T  K  STRM 
DIAGNOSTICS  OFF 
MONITOR  2  32  129  U 
MONITOR  4  32  1  U 
TTMESAVE  250 
UNIQUE_NAME  ON 

END 
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